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1 Summary 

Several numerical solutions of the three-dimensional unsteady incompressible Navier-Stokes equations in gener- 
alized coordinate systems are presented in this work. The governing equations are discretized by finite volumes 
with special care to the accurate approximation of the geometric quantities. The unknowns are the pressure and 
the volume fluxes over the computational cell faces. This formulation results in a robust fractional-step solution 
method for solving the discrete equations. Although the method is formulated for the three-dimensional case, only 
two-dimensional unsteady results are given in this work. Results are presented for a lid-driven two-dimensional 
cavity flow at Reynolds number of 10 4 , and for the flow over a circular cylinder with vortex shedding for several 
Reynolds numbers in the range 100 < Re < 1000. 

2 Introduction 

Numerous fluid flow phenomena are essentially time-dependent. Until recendy, the numerical simulation of unsteady 
flow fields over complicated configurations was very difficult because of limited computer resources both in terms 
of processor speed and storage. The new large supercomputers are powerful enough to make the unsteady flow 
simulation possible, provided an efficient and accurate solution procedure can be devised. 

This work presents a numerical simulation method for viscous time-dependent incompressible flow fields in ar- 
bitrary geometries. Two principal approaches have been used for the solution of the incompressible Navier-Stokes 
equations in primitive formulation [ 1 ] . One method is based on the fractional-step solution procedure while the other 
uses the artificial-compressibility concept with dual-time stepping for restoring the time-accuracy of the solution. 
The present study uses a fractional-step solution method. An efficient and robust solution procedure has been devel- 
oped and several test cases have been computed to validate the method. The details of the numerical procedure and 
some of the validation cases are presented in Refs. [2] and [3]. Therefore, only a brief description of the method is 
given here and a few solutions are included. 
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3 Formulation 


Navier-Stokes equations govern the flow of isothermal, constant-density, incompressible fluids. For a control volume 
V with surface S the mass conservation equation is 


i 


<2S • u = 0 , 


( 1 ) 


and the momentum conservation equation is 


4 [ udV=<f dS T, (2) 

at Jv Js 

where u is the velocity vector, t is the time, dS is a surface- area element, and dV is a volume-element. For Newtonian 
fluids the tensor T is given by 

T =-uu-Pj + jy(Vu + (Vu) r ) , (3) 

where P is the pressure, v is the effective kinematic viscosity, I is the identity tensor, Vu is the gradient of u, and 
(■) T is the transpose operator. 

4 Discretization and Numerical Solution 

Equations (1-3) are directly discretized over finite volumes to yield a fully conservative second-order scheme. Finite- 
volume discretization is preferred since it usually results in more accurate and stable solutions for generalized co- 
ordinate systems, especially for meshes with heavily clustered points and large curvature. Special care is given to 
satisfy the “geometrical conservation laws” [2] to minimize the errors resulting from the spatial discretization. 

The numerical problems associated with the absence of a pressure time-derivative in the mass conservation 
equation are handled by a fractional-step procedure. In each time step, the momentum equations are solved for 
an approximate velocity field which does not satisfy the continuity equation. In the second stage, the velocity and 
pressure fields are corrected such that the mass conservation equation is satisfied. This step leads to a Poisson 
equation with Neumann-type boundary conditions which may exhibit very poor convergence properties, especially 
in generalized coordinate systems. The CPU time consumed by the Poisson solver may be as high as 80% of the 
total computational time even in a Cartesian case [4], In the present work an attempt is made to minimize these 
difficulties. 

This attempt has led to the choice of the pressure (defined at the center of the computational cells) and the volume 
fluxes across each computational cell face as the dependent variables instead of the familiar Cartesian components of 
the velocity. (The volume fluxes are equivalent to the scaled contravariant velocity components in a staggered grid.) 
This choice ensures the satisfaction of the discrete mass conservation equation to round-off errors in any coordinate 
system and has favorable effect on the convergence properties of the Poisson solver (see details in Refs. [2] and [3]). 

A consistent Poisson solver is obtained by deriving the Laplacian operator from the discrete equivalent of the 
operator V • V (where V • and V are the divergence and the gradient operators, respectively), rather than discretizing 
the Laplacian operator directly. A novel and efficient ZEBRA scheme with four-color ordering is devised for the 
efficient solution of the nonorthogonal Poisson equation on vector computers. The Poisson solver could be converged 
to any specified small error in all the cases solved so far. 
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5 Results 


Various test cases were solved to validate the present solution method. The three-dimensional flow in a cubic cavity 
as well as the symmetric flow over a circular cylinder at Re = 40 and at Re = 200 with vortex shedding are described 
in Ref. [2]. The present paper reports additional cases to assess further the capabilities of the method. All the cases 
were solved using the three-dimensional code. Two-dimensional solutions were obtained by specifying along the 
third direction two intervals and periodic boundary conditions. 

5.1 Lid-driven Cavity, Re = 10 4 

In this case the time-dependent flow in a lid-driven square cavity was solved for a Reynolds number of 10 4 . A 
Cartesian nonuniform grid of 81x81x3 mesh points was used with clustering near the walls. The minimal interval 
is 0.28 • 10 -2 and the maximal interval at the center is 0.23 • 10 -1 . The upper lid started moving impulsively at T = 0 
and the time-accurate flow field was computed until the nondimensional time T = 120 was reached. Previously, 
steady solutions of this case were reported in the literature, see for example Ghia etal. [5]. However, these solutions 
were not obtained by a time-accurate method. For Re = 10 3 , the steady solution obtained by the present solution 
procedure compares favorably with other computations (see Ref. [3]). 

For the case Re = 10 4 the evolution of the main vortex and of the secondary vortices at the comers of the cavity 
exhibit complicated patterns that are repeated almost periodically without approaching a steady state. Figure 1 shows 
the instantaneous streamlines for 102 > T > 92 which corresponds approximately to one “period.” During this time 
interval, the main vortex rotates around a fictitious center and the flow field at the lower two comers goes through a 
complicated series of states which include the generation of secondary and tertiary vortices and their merging with 
one another. Recently, Bruneau and Jouron [6] have found that for this case the solution has not converged to a 
steady state if a fine grid is used. However, no experimental results are available for comparison at this time. 

5.2 Symmetrical Flow over a Circular Cylinder 

The symmetrical flow field over a circular cylinder at Re = 40 was solved employing a nonorthogonal grid with 
45 x 73 x 3 mesh points in the radial, circumferential, and axial directions, respectively. Figure 2 compares the 
time-evolution of the separation length with other numerical and experimental results. Good agreement is obtained, 
especially for T < 6. Slight disagreement is found for T > 6 due to the “wall effect,” which is the result of the 
finite distance of the outer boundary and the Dirichlet-type boundary conditions given in the present solution. This 
wall effect was studied earlier by varying the outer boundary distance (see Ref. [2]). 


5.3 Vortex Shedding from a Circular Cylinder 

This case considers the shedding of vortices behind a circular cylinder at Reynolds numbers 100-1000. A cylindrical 
grid of 81 x 84 x 3 mesh points was employed with points clustered near the cylinder and in the wake region. The 
vortex shedding was triggered by rotating the cylinder for a short period of time in both clockwise and counterclock- 
wise directions. Figure 3 shows the time evolution of the lift coefficient (dashed line) and the drag coefficient (full 
line) for several Reynolds numbers. In each case, the onset of the vortex-shedding and the periodic flowfield can 
be observed after a relatively short transient flow. Figure 4 compares the Strouhal number with other experimental 
and numerical results given in Refs. [9]-[ll]. The present solution yields a somewhat high Strouhal number for 
all the range of the Reynolds numbers. This behavior was also found by Gresho et al. [10]. Figure 5 shows the 
instantaneous streamlines for Reynolds numbers 100, 200, 550, and 1000 at the instant which corresponds to the 
maximal lift coefficient. 
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6 Concluding Remarks 


A time-accurate method for the solution of the incompressible Navier-Stokes equations is developed based on a 
fractional-step method. The dependent variables are the pressure and the volume fluxes across the computational 
cells. The governing equations are discretized by the finite-volume method with special care to satisfy the geometric 
conservation laws. The method is second-order-accurate both spatially and temporarily. 

The computed results compare favorably with other numerical and experimental studies. In the present paper 
some of the numerical results for the flow in a lid-driven square cavity and over a circular cylinder are presented. 
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Figure 2: Time-evolution of the separation-length (Re = 40) 
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Figure 3: lime-evolution of the force coefficients 
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Figure 4. Dependence of the Strouhal number on the Reynolds number (circular cylinder) 
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Figure 5: Effect of the Reynolds number on the instantaneous streamlines 
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